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CHEBYSHEV SERIES SOLUTION OF NONLINEAR ORDINARY 
DIFFERENTIAL EQUATIONS: INITIAL- VALUE 

PROBLEMS 

By Kin L. Lee and Paul F. Byrd 
Ames Research Center 


SUMMARY 


The approximate Chebyshev series solution of nonlinear ordinary 
differential equations based on Picard iteration is discussed. Detailed 
algorithms are provided for the numerical solution of initial-value problems 
involving (a) a system of n first-order nonlinear differential equations, 
and (b) an nth-order nonlinear differential equation. Methods to accelerate 
the convergence of the iterative procedures are also proposed. FORTRAN IV 
subroutines for the algorithms with options for accelerating convergence are 
given. 


INTRODUCTION 


Chebyshev series have been used for obtaining efficient numerical 
solution to various problems, particularly in approximating functions and in 
solving certain linear differential equations and Fredholm integral equations 
(refs, 1 and 2). The usefulness of the Chebyshev series lies in the fact that 
it very nearly satisfies the condition for optimal polynomial approximation. 

The explicit representation of the solution to a nonlinear ordinary 
differential equation by some economical approximation such as a Chebyshev 
series is particularly desirable in practice if the solution is part of a 
larger computation and is repeatedly required for the generation of other 
relevant information. Besides providing an explicit representation of the 
solution over the relevant range of the independent variable, the Chebyshev 
series solution of ordinary initial-value problems has other desirable fea- 
tures that may make it preferable to discrete variable methods (i . e . , Runge 
Kutta and various predictor- corrector procedures). First, the maximum error 
over the entire interval of integration can be readily and closely estimated 
by inspection of the coefficients. Second, since the integration can usually 
be effected over the entire interval or a small number of subintervals, the 
chances for propagation of round-off error is small. The degree M of the 
approximating polynomial as well as the interval of integration may be varied 
as subroutine parameters to obtain any degree of accuracy. These features are 
lacking in most discrete variable methods. A clear disadvantage of the algo- . 
rithms for Chebyshev series integration of nonlinear differential equations is 
the time-consuming feature of their interative construction. However, in the 
case of a single nth-order differential equation, the application of the 
Chebyshev series method in conjunction with the method of accelerating conver- 
gence may compare favorably with discrete variable methods in computing time. 



Approximate Chebyshev series solution for nonlinear differential 
equations was first proposed by Clenshaw and Norton (ref. 3) . The approximate 
solution, found by Picard iteration, is applied to single first and second- 
order differential equations. In a later paper (ref. 4), Norton proposed 
Chebyshev procedures based on Newton iteration to solve the same equations. 

The principal objective of this report is to give detailed algorithms 
based on Picard iteration for obtaining numerical solutions, in the form of an 
approximate Chebyshev series, of initial-value problems involving (a) a system 
of n first-order nonlinear differential equations, and (b) an nth-order non- 
linear differential equation. These algorithms have been programmed as 
FORTRAN IV subroutines. The documentation of the subroutines is provided in 
appendix C. 

A basic difference between algorithms presented here and those of 
Clenshaw and Norton is the choice of the interpolating polynomial. Our pref- 
erence is based on the discussion of error given in appendix B where certain 
basic ideas and tools of polynomial approximation are discussed. 

Although Chebyshev procedures based on Newton iteration usually seem to 
converge faster than those based on Picard, they are not easily extended as 
general algorithms to solve higher order or coupled differential equations. 

For this reason, no discussion on Newton's method is given here. To compen- 
sate for the slower convergence of Picard iteration, two methods are proposed 
in this report to accelerate the convergence of the above algorithms. The two 
methods are included as options in the subroutines mentioned earlier. 

A discussion on convergence via an example is given to provide insight 
into the behavior and application of the algorithms of this report. 

A summary is given in appendix A on several fundamental properties and 
tools of the Chebyshev polynomial. 


SYMBOLS 


Afc Chebyshev coefficients of the function <J>(t) (see eq, (Bll)) 

B k coefficients of the interpolating polynomial P M (t) (see eqs. 

(B23) and (B27)) W+1 

coefficients of the approximating polynomials for <j>4 (t) and <j>j(t) 
K (see eq. (14)) J J 

B^ 1 -^ coefficients of the approximating polynomials for 4. (t) and 


<f>v + 1 (t) (see eqs. (16) and (20)) 
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bf- 0 ?, bP? coefficients of the approximating polynomials for <j>. . (t) and 
X,K 1>K 

*(t) (see eq. (46)) 

1 9 J 

B? 0 ? , bP? coefficients of the approximating polynomials for <f>, . (t) and 

1 , K 1 , K X j jT'i 

<^j +1 (t) (see eq, (47)) 

bP^ coefficients of the approximating polynomial for <|>P^(t) (see eqs. 

K (67) and (68)) 3 

Bp^ coefficients of the approximating polynomial for $ CO (see eqs, 

K (67) and (68)) 3 

D n the set of polynomials of maximum degree n (see eq. (B2)) 


E N Cf) 

the minimax of f(x) (see eq. (B3)) 


P(x) 

a polynomial in x 


p k M 

a polynomial in x of degree k 


% CO 

the first N+l terms of the Chebyshev 
(B13) ) 

series of <j>(0 

T k Ct) 

the Chebyshev polynomial of the first 
(Al)) 

kind of degree 

“(it*) 

the maximum error of S N (t) (see eq. 

(B14)) 

^(t) 

the jth approximation to .(t) 


♦ P (t) 

the jth approximation to <{P-pt) 


A< 

E U k 

k=o 

= I U 0 + ui + . , . + U N 


n 




k=o 


2 U o + U 1 + 


+ U N^1 + 2 U N 


the convergence criterion or prescribed convergence error of 
algorithms I and II 

approximately equal to, as P M (t) « 4 (t) 
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PRELIMINARY ANALYSIS 


Chebyshev Series Integration of the First-Order 
Differential Equation 

The basic ideas and tools used in the construction of algorithms for the 
approximate Chebyshev series integration of nonlinear differential equations 
can best be discussed and understood via the first-order differential equation 

§=fCx,F) (1) 


having the initial condition 


F (a) - n (2) 

(Algorithms for a system of first-order differential equations and an nth- 
order differential equation are presented in the following sections as algo- 
rithms I and II in a form suitable for coding by means of ALGOL, FORTRAN or 
similar computer languages.) 

Consider the sequence of functions (Fj(x)} generated by a process 
attributed to Picard: 


F, + 1 (x) = n + J X f (s,Fj)ds (j = 0,1,...) 


(3) 


FqOO = T) (4) 

If f(x,F) is continuous and the partial derivative 8f/£F is bounded in a 
region including the point (a,n), the above sequence of functions is guaran- 
teed by a theorem of Picard to converge to a function F (x) satisfying 
equations (1) and (2) in a neighborhood | x - a | < h. If this is the case, 
then without loss of generality we can consider the solution of the differen- 
tial equation 


ft = » 


-1 < t < 1 


(5) 


with 


<K-l) = n 

by means of the iterative procedure involving the equations 


ft 

<f>j +1 (t) = n + J iKu,<f>j)du 


( 6 ) 


(7) 
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3 n (8) 

where ^ (t) converges uniformly on [-1,1] to the solution of equations (5) 
and (6) . 

Although i(»(t,<(ij), for a fixed j, is an explicit function of t, the 
integral of (7) is difficult to obtain in practice. However, if (t)] 

can be accurately approximated by a polynomial P(t), then <l>j + iCt) can be 

evaluated by integrating P(t) term by term. From the point of view of effi- 
cient computation, the coefficients of such a polynomial should by readily 
obtainable by a finite algorithm. Also, for a fixed degree M, this polynomial 
should be the best possible in the sense of least maximum error (defined by 
eq. (B2)). Since the computation of the best approximating polynomial is, in 
general, a nonlinear iterative procedure, the use of it as an effective tool - 
in the approximation of iHt,<f>j) must be precluded. Clenshaw (ref. 3) pro- 
posed the use of the interpolating polynomial 1 

JLh 

PmOO = T, C k T k (t) (9) 

k=o 

where T k (t) are Chebyshev polynomials defined by (see also appendix A) 

T k (t) = cos(k cos -1 t), -1 < t < 1 (10) 

with the points of interpolation 

TTT 

t r = COS S. (r =0,1, . . . ,M) (11) 

where T M (t) has M + 1 extrema T M (t r ) = (-l) r . 

Here, however, we make use of Qjq(t) , a modified interpolating polynomial, 
which is formed by truncating the last term of the interpolating polynomial 

M+l 

W« ’ S" C 12 ) 

k=o 

having 

t r = cos gL_ (r = 0,1, . . . ,M +1) (13) 

as the points of interpolation. These are also the M+2 points where 

Tm+i ( t) has extrema, % +1 (t r ) = (-l) r . Analysis in appendix B shows that the 

maximum error for %(t) as an approximating polynomial for sufficiently large 

~ ^A double prime over the summation sign - indicates that the first and last 
terms are to be halved, while a single prime indicates that only the first 
term is to be halved. 
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M is one half that of P^ft). Numerical examples of Frazer and Hart (ref, 5) 
also show that Q M (t) closely approximate the best approximating polynomial. 
For this reason, Q^(t) is also called a near-best approximating polynomial. 

Now suppose we assume that each member of the sequences (<J>j(t)} and 
U j (t) } can be accurately represented by polynomials of degrees M + 1 and M, 
respectively. Assume also that at the jth iteration <j>j (t) and <j>j (t) are 
known and of the form 


M+l 


1 


= E' b k 0,T kW 

k=o 


and 


(14) 


M 


♦jco b k l)T kW 
k=o 

'then an approximate chebyshev series solution for equation (5) can be obtained 
as follows: 

We approximate first <f>j+i(t) = d* [t , 4> j (t) ] by the near-best approximating 
polynomial Qm ( t) to obtain 

M 

HtAj (t)] -X) B k° T k(t) (15) 

k=o 

where according to equation (B27) is 

M+l 


B k ^ M+l ^ ^t t r»^j ^r^-r^k^ 
r=0 


(16) 


with 


„ ru 

t_ = cos rr— 
r M+l 


(17) 


Since B^ 1 ^ is a linear combination of the form (A22), it can readily be 
evaluated for a fixed k by means of recurrence formula (A23) , Accordingly, 

C M+1 = 2 ^^M+l ^M+l-^J 

C M " ^f t M , < * > j ^M-O + 2t k c M+l 

c r = (*r)J + 2t k c r+l " c r+2 

(r = M - 1,M - 2, v > .,1) 


(18) 
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4 ^ = I ,f) [ t o’ ()> j (t o^] + ®l l k " c 2 


(19) 


Denote the integral of equation (16) by the (M+l)st degree polynomial 

f M , f N M+l f N 

* j+ ict) - fx: 1 - X’ 4 

J i „ V-n 


Upon performing the indicated integration with the aid of equation (A9) and 
equating coefficients of T^(t), one obtains 


R Cl) B (l) 

■ Co) _ k-l ~ Vi 


(k = 1,2, . . . ,M) 


a (0) _ d m 

M+l M+l 


(The same results can also be obtained if equation (B19) is applied.) The 
constant of integration (1/2) B^ remains to be determined. It can be 


computed if one notes that 


and by equation (A31) that 


<t>j + l(-l) = h 






B o o) ■ 2 [ n -2^c-i) k ^ o) ] 


This completes one iteration. If 


b k l} " 4^ < e Ck = °» l » ‘ ' -’ M) 


where e is a prescribed convergence error, we are through. Otherwise, 
replace each b£^ by B^ (p = 0,1) and initiate another solution.. 
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The entire iterative process can be started by setting <f> 0 (t) = n» that 
is, by taking b^°^ = 2n and b/°^ = 0, for k = 1,2, . . . ,M+1. A better 

O K 

initial approximation sometimes can be made by utilizing knowledge of ik(t,<|>), 
for example, if ip (t,<j>) involves only the dependent variable <f> . 


Accuracy of Solution 


When the condition (24) is met, we are interested in how well the 
approximate solution satisfies the given differential equation. If the poly- 
nomial approximations (15) are substituted into equation (5), one obtains, 
upon taking absolute values, 


M 


b k\W 


- t 


k=0 


M+l 


*E 


L k=0 



M M 


M 

_ „ f 

~ M+l 


< 

E' b k° T k (t) -E B k 1)T k(o 

+ ■: 

E »k 1)T kto - * 

c E' b k\v 



| k= o k=0 


k=o 

- k=o 



But 


1 k=o 

by condition (24) . Let 


4> 


M 

; b k 1 )T kt« -E B k l)T ict« 

k=0 


< (M + l)e 


M+l. 

E 


b ^ J T k (t) 


00 - 

E 


L k=o 

It then follows from equation (B36) that 


k=0 


M , f M+l -] 




00 , 

2 E^TfcCt) - ♦ t-E b k 0)T kW 


< 

41 ! 

•■EK" 

k=0 L k=0 




k=M+2 


(25) 


( 26 ) 
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Consequently, we see by equations (25) and (26) that 


b£ 0) T k (t) - i, 


< cm* vKid 


k=M+2 


Hence if b^ 1 ^ - B^ 1 I < e for each k and if M is sufficiently large, the 

M+l t ^1, 

approximate solution ^ ^ bj^Tktt), and, a fortiori, B k °^T k (t), will 

k=o k= o 

satisfy the given differential equation with an error close to (M + l)e. In 
practice, we say that M is "sufficiently large" when larger values provide 
no change greater than e in the coefficients. Also, since the Chebyshev 
series is unique, when (M + l)e is made small by an appropriate choice of e 

and M, the coefficients B^ 0 ^ of the finite series will closely approximate 
those of the Chebyshev coefficients of the solution (see eq. (Bll). 


Example 

To illustrate the accuracy of the above procedure, let us find a 
polynomial approximation for tan[ (*rr/8) (t + 1)] for -1 < t < 1 with a maxi- 
mum error less than 0.5xl0“ 8 . One can easily verify that the given function 
satisfies the differential equation 

gf- = f[l + -1 < t < 1 (28) 

with 

<K-1) - 0 (29) 

Hence, the method given in this section is applicable. 

The approximate Chebyshev coefficients for both the solution and its 
derivative corresponding to M = 16 and e = 0.5xl0 -1 ° are shown in table 1(a). 
A total of 13 iterations were required. Tabulated values of the approximate • 
Chebyshev series corresponding to discrete points of the independent variable 
are given in table 1(b). Numerical results suggest that M == 16 is suffi- 
ciently large since the coefficients corresponding to e = 0.5x10" 10 for 
M > 16 yield no change greater than e. In view of equation (26), the approxi- 
mate solution must satisfy the differential equation with an error bound close 
to 8.5xl0 -1 °. The same conclusion can be drawn by the examination of the coef- 
ficients alone. In fact, since b/ 1 ^ is approximately equal to a/ 1 ^ 
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< 1/5 for K > 7, the right 


member of inequality ( 27 ) gives us 

oo 

17e + ax 7 + Ja^l « 8 , 7x10 " 10 

k=i8 

A check of the tabulated values of table 1(b) with those of Abramowitz and 
Stegun (ref. 6 ) shows agreement to 10 decimal places. 

TABLE I .- CHEBYSHEV SERIES APPROXIMATION OF TAN j|(t + 1)1, (-1 < t < 1) 


(a) Approximate Chebyshev coefficients 


k 

d(o) 

B k 

R (l) 

B k 

0 

0.9113043408269388D 00 

0 . 1Q43307398148425D 01 

1 

0 . 4894686436450291D 00 

0. 1835797842596252D 00 

2 

0 . 4284834890908355D-01 

0. 6437011085836632D-01 

3 

0. 1024434335477160D-01 

0. 1218638862329105D-01 

4 

0 . 1453268753787471D-02 

0 . 2904050729736723D-02 

5 

0 . 2787802534394446D-03 

0. 56023859299 12773D-03 

6 

0 . 4483018228371541D-04 

0. 1162481953422768D-03 

7 

0. 7992572697332448D-05 

0 . 2227640558669249D-04 

8 

0.1340847727289307D-05 

0.4352177579622525D-05 

9 

0 . 2331260728730840D-06 

0.8228419500635855D-06 

10 

0 . 3968754154824808D-07 

0 . 15590826790701 23D-06 

11 

0 . 684067079066359 ID-08 

0. 2909111909862381D-07 

12 

0 . 1170504525821132D-08 

0.5413510512413251D-08 

13 

0 . 2011467240884598D-09 

0. 9990104789166475D-09 

14 

0 . 3447781322812264D-10 

0. 1836956861132969D-09 

15 

0 . 5912174018325965D-11 

0 .336317085292 1371D- 10 

16 

0. 1050990891537928D-11 

0 . 6330465563517943D-11 

17 

0. 1861901636328807D-12 

0 . 0000000000000000D-38 


according to equation (B28) and 


(b) Function values 


t 

tan[j-(t + 1 )] 


- 1.0 

O.OOOOOOOOOOOOOOOOD- 

■38 

- 0.8 

0 . 7870170682457329D- 

-01 

- 0.6 

0 . 1583844403245379D 

00 

-0.4 

0 . 2400787590801460D 

00 

- 0.2 

0 . 3249196962328421D 

00 

0.0 

0.4142135623731530D 

00 

0.2 

0 . 5095254494943512D 

00 

0.4 

0 . 6128007881399821D 

00 

0.6 

0 . 7265425280053249D 

00 

0.8 

0 . 8540806854634090D 

00 

1.0 

0 . 9999999999998522D 00 
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For economy of computation, note that if Ni < N, then 


N t , 

2J B k T k (t) - 2 B k T k (t) 
k=o k=o 


(See also eq. (B16)). Thus, in the case of the approximating polynomial of 
the above example, ignoring the last six terms from the finite series results 
in a maximum error of 0.14xl0~ 8 . Hence, an 11th degree instead of a 17th 
degree, polynominal can be used to approximate tan[ (tt/3) (1 + t)] and still 
satisfy the maximum error requirement of O.Sxl.0' -8 . 



INTEGRATION OF A SYSTEM OF n FIRST-ORDER DIFFERENTIAL EQUATIONS 


In this section, the basic ideas applied to the construction of an 
approximate Chebyshev series solution for a single first-order differential 
equation are extended to provide an algorithm for the solution of a system of 
n first-order differential equations. It is given in sufficient detail to 
facilitate computer programming as well as. the discussion of acceleration of 
convergence. The basis for the more general algorithm is the following 
theorem . 

•* 

Theorem . Let a system of n first-order differential equations be 
defined by 

= fjCx,Fi,F2> • « • » Fj^) , jx - a | < C Q , (i — 1,2, . , . ,n) ■ (31) 

with the initial conditions 

F t (a) = n ± , Ci - 1,2, , . . ,n) (32) 


Furthermore, let each of the functions f. be continuous and have bounded 
partial derivatives 


8f 

9F 


1 


< K , (i, j = 1,2, . . . ,n) 


(33) 


in the region 


x - a < 


< C 0 , F^ — rii - , (i — 1,2, . . . ,n) 


(34) 


If 


h 


min^C 0 , 


Ci 

L 


C 2 

T 



(35) 
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where 


L = max£nax|f L (x,F 1 ,F 2 , . . . ,F n )jj 

in the region defined by (34) , then the sequence of n functions 

(pij W,F 2> j(x), . . .,F n ,(x)| defined by 
' lj=0 


(36) 


F i, j + i « = + J a f i^ s ’ F i, j >^ 2 ^ j » • • (i = 1»2, . . . ,n) 


(37) 


with 

F i>0 (x) = n 1 , (i = 1,2, . . . ,n) (38) 

converges uniformily on |x - a| < h to a unique set of functions of Fi(x), 
F 2 (x), , , F n (x) satisfying equations (31) and (32). (For proof of a simi- 
lar theorem, see Tenenbaum and Pollard (ref. 7)). 


Besides providing an iterative procedure (eqs. (37) and (38)) for 
obtaining a solution, the above theorem also guarantees an interval of conver- 
gence. However, the estimate h in equation (35) usually proves to be con- 
servative if not difficult to find. In practice, the interval of convergence 
is usually assumed or determined by trial and error. 


Suppose that the system (31) has a unique solution F^(x) 

(i = 1,2, . . . ,n) on a < x < b satisfying the initial conditions given by 
equation (32). In order to construct the Chebyshev series for F^(x), make 
the change of independent variable 

x = ct + d, 

with 


and 


c 


b-a 

2 




(39) 


so that 


b+a 

2 


J 


F t (x) = <f>i(t) , 



c"Hi(t) , 


F t (a) = ^(-1) = m 


Substitution in equations (31) and (32) then yields 

Ct) = Ct, 4>i ,4>2, • • • <t> n ) > -1 < t < 1 , (i = 1,2, . , . ,n) (40) 

with 


<i>i(-l) = ni , (i = 1,2, . . . ,n) 


(41) 
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The same change of variable for equations (37) and (38) gives the sequence of 
n functions |(j>i > • • * » 4*n , j } defined by 


4>i, j+1 Ct) hi '{'i ( u > ^ i , j » >4> n , j)du j (i = • « • >n) 

_1 (42) 

and 

<t>i }0 (t) = n ± » (i = 1,2, . . . ,n) (43) 

which converges uniformly on the closed interval [-1,1] to a set of functions 
<h(t)> ^(t), • • •» <t> n (t) satisfying equations (40 and (41). One also 
obtains by differentiation the equations 


with 

= ^(t,ni,H2, 


(i = 1,2, . . . ,n) 
.,n n ) 


(44) 

(45) 


We are now ready to proceed with the algorithm for an approximate 
Chebyshev series solution of n first-order differential equations. 


Algorithm I 

As in the case of a single first-order differential equation, it is 

assumed that the sequences!^ and 

■ j = 0 

j't’i j>'$2 j > • • •j'f'n )| can be accurately approximated by polynomials of 
degree M+l and M, respectively. For simplicity of notation, let 


M+l , 

Z b k!i T k (t) ’ 

M , 

b M T k tt5 

(46) 

k=o 

k=0 


M+l 

M 


Z' - 

*k!Kw 

(47) 


k=0 


k=o 
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Also, let e be a prescribed convergence error between and B^ 1 }. The 

approximate solution of (40) can be obtained as follows: 

1 • Set kTr 

= cos » k = 0, 1, . . . , M + 1 
(These are the M + 2 points where T m+ ^Ct) assumes its extrema.) 

2 . Set 

(n'\ 

b^ 


= 2n i , 

(i = 1,2, . . 

• >n) 

= o , 

(i = 1,2, . . 

• • • • jM + 

= 0 , 

(i « 1,2, . . 

• • o • 


k,i 

b^ 1 ) 

k,i 

(This is equivalent to the initial approximations d>i (t) = n^, 
(i = 1,2, . . . ,n) ) . 

3 . Compute 

M+l 




(i = 1,2, . . ,,n; k= 0,1, . . . ,M + 1) 


r=0 


4 . Compute 

(i = 1,2, , . . ,n; k = 0,1, . . . ,M + 1) 

5. Compute (by making use of eq. (B27)) 

M+l 

B k!i = M+T ^ *i,j + i^WV * (l = 1 > 2, * * * ,n; k=0,1, • . . ,M) 

0 r=o 

(We have here approximated <(>! , (t) by the Mth degree polynomial 

i , 1 + 1 

M 


B^^Tkft) , where according to Theorem B5 


k=0 
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B k!i • • * ^n,j] T kCt) Cl - 1 2 r 1/2 dt 

for a sufficiently large M.) 

6. Compute (by eq. (B19)) 

. _ jjCO 

B ( 0 ) = k ~ 1 » i 2K k+1>i , (i = 1,2, . . . ,n; k- 1,2, . . . ,M) 


M+l , i 


1+11 * ^ ~ ‘ 


7. Compute (using the fact that <j>. . (-1) = n- and eq. (A31)) 

1 9 3 + 1 1 


B*- 0 } = 2 
o,x 


M+l 1 

-i ■ 


(i - 1,2, • * • >n) 


(We have in this and the above step integrated 


i v .i 

i.M w 


which results in the (M + l)st degree polynomial 

M+l , 

- X) 


If IB^ 1 } - b^ 1 ?! < e for all k and i, we are through. Otherwise, 
K ) 1 K , 1 


9. Set 


b k?i * b £1 • tp - 0,1) 


for each k and i and return to step 3. 

Note: To take into consideration the case when <j>! (t) = 

i, step 8 should be bypassed until the second iteration. 


1>! n (t) = 0 for each 
1 , u 
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Example I 


As an example, we solve the boundary- layer equation 


d 3 F 

dx 3 


+ 2F 


d 2 F ( dF 


-(g) *1 = 0 , 0<*<6 


( 48 ) 


with the initial conditions 


FCO) = 0, gj 


x=o 


= 0 , 


d 2 F 

dx 2 


= 1.311937693880 


(49) 


x=o 


using algorithm I with a prescribed convergence error e = 0.5x10” 10 . 

We first rewrite the given differential equation as a system of three 
first-order differential equations. Let Fi(X) = F(x), F 2 (x) = dF/dx, and 
F 3 (x) = d 2 F/dx 2 . Equation (48) can then be equivalently written as 


dF x dF, dF. 

Ir = F 2 (x), 2 ^= F 3 (x), ^ = F 2 2 (x) - 2F 1 (x)F 3 (x) - 1 

with initial conditions Fj{0) = 0, F 2 (0) - 0, F 3 CO) = 1.311937693880. 

We were unable to find a solution to the example by the method of this 
section over the entire interval because the iteration process failed to con- 
verge. However, we can obtain a solution over subintervals of the given 
interval. Let the given interval be subdivided by the q + 1 points 

0 = Xp < xj < . . , < Xp = 6 (50) 

On each of the subintervals [Xy.^Xj.] (r = 1,2, . . . ,q) the change of 
variable x = ct + d, c = (x r - x r _ 1 )/ 2 , d = (x r + x r _ 1 )/2 enables one to 
rewrite equation (48) as 

4>1 (t) = c<t > 2 (t) 2 , <j> 2 (t) * ctjj^Ct) , <t> 3 (t) = c[<j> 2 2 (t) - 2 ij)j (t)<t> 3 (t) - l] 

The solution can then be found for one subinterval at a time. Initial condi- 
tions for each subinterval are the function values of the end points of the 
previous subinterval, except the first where 4>i (-1) = Fi(0), 4 > 2 C-l) = F 2 (0), 

and <j> 3 (-l) = F 3 (0) . The approximate Chebyshev coefficients of 

Fi(x) = F(x) for the case with M = 11 and x Q = 0, xi = 1, x 2 = 2, . . . , 

Xg = 6 are given by table 11 (a). The coefficients for the first and second 

derivatives are not tabulated but may be generated in terms of from 

K, 1 

equation (B19) . A tabulation of the values of F(x), dF/dx, and d 2 F/dx 2 
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corresponding to x •= 0(0. 2)6.0 is provided by table 11(b). Numerical results 
indicate that the approximate solution satisfies the differential equation 

with an error bound of the order O.SxlO" 1 ^, since exhibits no change to 

K,X 

10 decimal places for M > 11 (see Preliminary Analysis). Note that the first 
and second derivatives approach unity and zero (accurate to 10 decimal places) , 
respectively. This is because the initial value of the second derivative was 
determined numerically from the solution of the boundary-value problem 
involving the same differential equation (45) and the boundary values 


F(0) = 0 , 


dF 

dx 



dF 

dx 


X=oo 


1 


(See L. Fox (ref. 8) for the numerical solution of boundary-value problems by 
means of initial-value techniques) . 
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TABLE II.- APPROX IMATATE CHEBYSHEV SERIES OF EXAMPLE I; 
M = 11 AND e = 0.5x10' 10 


(a) Approximate Chebyshev 


0. 3892366850285570D 00 

1 0 . 2510773451273367D 00 

2 0. 5 14969731 7862475D- 01 

3 -0.4884131498391032D-02 

4 0.9 1475 38219 192970D-04 

5 0. 1415676089930017D-04 

6 0. 4392922674440762D-06 

7 -0. 1455950532358585D-06 0 < x < 1 

8 -0 . 4849170368489048D-08 

9 0. 8282513783615847D-09 

10 0 . 9882967257068829D-10 

11 -0. 6903978861669351D-11 

12 -0. 8827151972331154D-12 

0 0. 1906417059176552D 01 

1 0. 472 151 77725 131 79D 00 

2 0.9310264019825763D-02 

3 -0. 1833951624942309D-02 

4 0. 2059970855772206D-03 

5 -0 . 8457783925547068D-05 

6 -0. 8661141518947492D-06 

7 0. 1097311409086846D-06 1 < x < 2 

8 0 . 3074443225322085D-08 

9 -0 . 1144396357769898D-08 

10 0. 1743676140975036D-10 

11 0.9458996086397773D-11 

12 -0.4752447471050794D-12 

0 0 . 3863311229682918D 01 

1 0. 4991658557862458D 00 

2 0. 3890170133408291D-03 

3 -0 . 1226429401112712D-03 

4 0. 2723834050810958D-04 

5 -0. 4276208831696196D-05 

6 0 . 4483551620850419D-06 

7 -0. 2359560210099267D-07 2 < x < 3 

8 -0 . 1071584776352978D-08 

9 0 . 3086935209765476D-09 

10 -0. 1991941542267561D-10 

11 -0. 7958215315696075D-12 

12 0 . 2178984230703938D-12 


coefficients of F(x) = Fi(x) 


IV, i 

0 0 . 5862202712178695D 01 

1 0 . 499994851464839 7D 00 

2 0. 2955377779468674D-05 

3 -0. 1245878270158787D-05 

4 0.40055 165615 12526D- 06 

5 -0. 1007087050914716D-06 

6 0. 200679 101 7080924D-07 

7 -0. 3171829281200812D-08 3 < x < 4 

8 0. 3907806881905951D-09 

9 -0. 3545559422 18950 ID- 10 

10 0. 1928744432685055D-11 

11 0. 2469520799973746D-13 

12 -0. 1720633666410817D-13 

0 0 . 7862196443651593D 01 

1 0 . 4999999943549484D 00 

2 0. 3692652891795766D-08 

3 -0. 1878812574654078D-08 

4 0 . 7638849269566333D-09 

5 -0. 2536970874119266D-09 

6 0. 7001695230546277D-10 

7 -0. 1625416567214393D-10 4 S x s 5 

8 0 .3200698242 324569D- 11 

9 -0.5361639630217125D-12 

10 0. 7642312708592647D-13 

11 -0.9226205658049735D-14 

12 0.9246076126956381D-15 

0 0.9862 196437 12 1186D 01 

1 0 . 5000000000008645D 00 

2 0. 7747691377346653D-12 

3 -0. 4426482328826845D-12 

4 0.1981233464324698D-12 

5 -0. 8628422050923251D-13 

6 0. 2645800245559826D-13 

7 -0. 9037083251041076D-14 5 < x < 6 

8 0 . 1891137709394048D-14 

9 -0. 572 844241 17811 88D- 15 

10 0. 7540264708912517D-16 

11 -0 . 2460153293203471D-16 

12 0.2698458740408366D-17 
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TABLE II.- APPROXIMATE CHEBYSHEV SERIES OF EXAMPLE I; 
M = 11 AND e = 0.5x10“ 10 


(b) Values of the approximate solution and the first two derivatives 

for x = 0(0. 2)6.0 


\ 

X 

F (x) 


dF 

dx 


d 2 F 

0.0 

0.1387778780781446D- 

16 

0.000000Q00O000000D- 

•38 

0. 1311937693880000D 01 

0.2 

0. 249056486064939 1D- 

■01 

0. 2423943538939113D 

00 

0.1112107122641016D 01 

0.4 

0.9430258133873859D- 

■01 

0. 44498662 13741513D 

00 

0.9145465290677766D 00 

0.6 

0. 2003061101255821D 

00 

0. 6087099439315713D 

00 

0 . 7245087590634542D 00 

0.8 

0. 3353391239337122D 

00 

0. 7357728939772173D 

00 

0. 5492 147378 176965D 00 

1.0 

0.49241445 1232278 ID 

00 

0 . 82986805 10454096D 

00 

0 . 3959361369563590D 00 

1.2 

0.6654189143037998D 

00 

0 . 8959772698742670D 

00 

0. 26999030965 75698D 00 

1.4 

0. 8493213259061367D 

00 

0.939 82673439779 7 ID 

00 

0 . 1733584569022934D 00 

1.6 

0.1040250891316518D 

01 

0. 9671738223714611D 

00 

0. 1044301711933930D 00 

1.8 

0.1235435703128323D 

01 

0. 983158 1570755 136D 

00 

0 . 5884905029559881D-01 

2.0 

0. 1433033404109585D 

01 

0 . 9918921621644430D 

00 

0 . 3095379362802304D-01 

2.2 

0. 1631909451187109D 

01 

0.9963447565741975D 

00 

0. 1517036399321181D-01 

2.4 

0 . 1831417188863473D 

01 

0 . 9984593547253356D 

00- 

0.6918247270753256D-02 

2.6 

0. 20312156729 18 176D 

01 

0.9993937528203719D 

00 

0 . 2932544923138562D-02 

2.8 

0. 2231138667864316D 

01 

0. 9997775510127270D 

00 

0.1154412461402519D-02 

3.0 

0. 2431111230808782D 

01 

0.9999239691430687D 

00 

0.42 172645422 12354D-03 

3.2 

0. 2631102124524431D 

01 

0.9999758156979835D 

00 

0. 1428864892765967D-03 

3.4 

0. 2831099311572934D 

01 

0. 9999928465 126406D 

00 

0 . 4487654564845384D-04 

3.6 

0. 3031098503447131D 

01 

0 . 9999980337347978D 

00 

0. 1305947442842973D-04 

3.8 

0. 3231098287668217D 

01 

0.9999994980713486D 

00 

0.3520034 178472342D-05 

4.0 

0. 3431098234149990D 

01 

0.99999988 10724133D 

00 

0. 8785038983643211D-06 

4.2 

0. 3631098221826554D 

01 

0.9999999738584204D 

00 

0 . 2029509793444975D-06 

4.4 

0 . 3831098219193559D 

01 

0 . 9999999946730537D 

00 

0.4338982542142417D-07 

4.6 

0. 40310982 186721 20D 

01 

0 . 999 999 99 89958 59 6D 

00 

0. 8583529424554144D-08 

4.8 

0. 4231098218576820D 

01 

0.9999999998268 14 ID 

00 

0 . 1572242307020353D-08 

5.0 

0.4431098218561269D 

01 

0.9999999999737648D 

00 

0. 2766056552012288D-09 

5.2 

0 . 4631098218559299D 

01 

0 . 9999999999991659D 

00 

0. 4395167933392727D-10 

5.4 

0 . 4831098218559648D 

01 

0. 1000000000003167D 

01 

0. 7029429704709523D-11 

5.6 

0. 5031098218560366D 

01 

0. 1000000000003867D 

01 

0. 1593260600584493D-11 

5.8 

0. 5231098218561158D 

01 

0 . 1000000000004142D 

01 

0. 3583023509204903D-12 

6.0 

0 . 5431098218561917D 

01 

0 . 1000000000005115D 

01 

-0 . 8557528262266267D-11 
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Remarks on Convergence 


In the application of the method of the foregoing section the integration 
involved in each of the iterations 







• ^n,j) du 


is an approximate one. The accuracy by which 'f’j., j + i (t) can be evaluated 
depends on how accurately ^(t,^ j ,<f> 2 , j , • • • j^n/i^ I s approximated, or, 


equivalently, on the degree M of the polynomial 



B^TkCt) 


used to 


approximate j , . , . , cfs^ , j !) * It is to be expected that the rate 

of convergence of the functions (t) may be influenced by M. Table III 

suggests that this is indeed the case. Numerical results showed that an 
M > 11 is sufficient to provide an approximate solution to 10 decimal places 
for example I, However, the number of iterations required, for a fixed conver- 
gence criterion e, may vary with M. Note that the number of iterations 
needed for the first two intervals remained fixed for each M. On the other 
hand, the number of iterations required for convergence over the remaining 
intervals decreases with increasing M, The same table shows, in the case of 
5 < x < 6, that there is a rather sharp decrease in the number of iterations 
when M is incremented only by 1. 


TABLE III.- THE INFLUENCE OF M ON THE CONVERGENCE OF ALGORITHM I. 


^Interval 

M 

0 < x < 1 

1 < x < 2 

2 < x < 3 

3 < x < 4 

4 < x < 5 

5 < x < 6 

11 

17 

21 

30 

38 

44 

53 

12 

17 

21 

28 

36 

41 

42 

13 

17 

21 

26 

34 

39 

39 

14 

17 

21 

25 

32 

37 

29 

15 

17 

21 

25 

29 

33 ! 

28 


Each entry of this table indicates the number of iterations required by 
algorithm I to solve equation (48") having the indicated M and interval of 
integration. The convergence criterion is e = 0.5xl0 -10 . Since the approxi- 
mate Chebyshev coefficients agree to 10 decimal places for all M > 11, the 
solution for each M is also accurate to 10 decimal places. 

An observation one can make is that with an appropriate choice of M the 
computing time used for solving a particular problem can be minimized. When 
an appropriate choice of M is not available it may be prudent to use a 
larger M than is deemed necessary for a prescribed accuracy (e.g. , using 
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M = 15 instead of 11 for example I) to insure convergence within a moderate 
number of iterations. In the next section it will be shown how the conver- 
gence of the iterative process can be accelerated for a fixed M, thereby- 
improving the efficiency of algorithm I, 


Acceleration of Convergence 


To economize computation, two methods consisting of "modification of rows" 
and "modification of columns" are proposed in this section to accelerate the 
convergence of the algorithm. Before we proceed further, the basic tool 
required by these methods must be given. 


Steffensen's sequences .- Consider first the modified sequence associated 
with the single-point scalar iteration function yj +1 = f(yp having the prop- 
erties: (a) the equation y = f(y) has a solution y = u; (b) the third deriva- 
tive f" ' (y) is continuous in a neighborhood of u with f"'(u) 1. Let the 

sequence be denoted by 



(51) 


The first member y^ 0 ^ of the sequence is an initial approximation to u. 
The other members are evaluated in the order indicated by the formulas 


yi (r) - f(/o (r) ) , y| r) » f(yi (r) ) 


(52) 


yo 


o (r+l) « < 


(y 2 (r) - yf r3 ) 2 


y 2 (r) - ^ - 1 - / , if y| r) - 2y 1 (r) + y 0 (r) + 0 

y| r) - 2y^ + y 0 Cr) 


(53) 


V W 

Yz 


if - 2,S r > * yo (r) = 0 


It can be shown that the sub-sequence |y having any prescribed yo is 

quadratically convergent (see, e.g, , ref. 9). Therefore, the application of 
equations (52) and (53) would effectively result in accelerating the conver- 
gence of the iteration function yj +1 = f(yj) whose regular sequence may be 

only linearly convergent. Equation (53) is known as Aitken's formula, but the 
scheme consisting of equations (52) and (53) in the evaluation of the sequence 
(52) is due to Steffensen (see ref. 9) . 
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^ Consider next the more general case of the single-point vector iteration 
f defined by 


yj + i = f( yp 

where y. is the N-dimensional vector 


(54) 


3 




■W 


it is also assumed that the vector equation y = f (y) has a solution y = u, 
and f(y) has certain desirable properties analogous to the scalar function 
f(y). However, even without knowing precisely what the desirable properties 
should be, we can, at least formally, generate a sequence of vectors analogous 
to the sequence associated with the scalar iteration function. This is accom- 
plished by recomputing a new after every N + 1 sucessive applications of 

(54). 


: be denoted by 


-* +(o) 

y o = yo 

+ (l) 
yo 

?o (2) 

Ho) 

Yl 

fy ( f‘ 

tvP 


1 rP 

^ 1 

/ 

1 

y N+l 

1 

yCO 

y N+l 

« 

H2) 

y N+l 


(55) 


Define two NxN matrices AYq and A 2 Yq by 
4Y 0 = 

a 2 y 0 » (t 2 yo r ^.4 2 yf r:) , . . . .4 2 y^i) 

with 

Ayf r) = y« - yW , A 2 ?? 5 = yf r > . 2yf r) * y« 

'j + l •’j * J+2 y j + l y j 

The first member y 0 of the sequence (55) is an initial approximation of u. 
The other members are evaluated in the order indicated by the vector 


22 



equations 


= f(yj r) ) (j = 0 , 1 , . . . ,n) (56) 

if det ^A 2 Yq^ ^ 0 

(57) 

if det ^A 2 Y - 0 

In particular, when N = 1, equations (56) and (57) reduce to (52) and 
(53) , respectively. 

Steffensen's iteration procedure for the vector case has not been 
investigated fully from the theoretical point of view. However, in practice, 

the sub-sequence jyo^j , for a large number of cases has been found, as in the 

scalar case, to be quadratically convergent (see ref. 9). Consequently, 
unless the regular sequence generated by equation (56) is already quadrat- 
ically convergent, it would be less rapidly convergent than that of 
Steffensen's. Even in the case when the former is divergent, the latter has 
been found to be convergent in many cases . 

The remainder of this section discusses the utilization of Steffensen's 
sequence in the acceleration of convergence for algorithm I. The approach 
here is equally applicable to algorithm II. 


yo Cr+1) =J 


►(r) 

N 


AY, 


i ( a 2 y o) 


- 1 A iCr) 


Ay. 


N 


, y W 

v N+ 1 * 


Recall that we have from algorithm I the sequences f(j>! . (t, )} 

\ K I j=o 

(i - 1,2, . . . ,n;k =0,1, . . . ,M+2), where each iterate of each sequence 
is generated at step 4 of the same algorithm. These sequences can be put in 
the form of a single sequence of rectangular arrays 


*1,1, j 


• • • 




• • • 

^2,M+2, j 

_fn,i ,j 

*n,2 ,j 

• ♦ • 

^n,M+2 , j 


(j = 0,1,2, . . .) 


(58) 
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where 




The (M + 2) column, being <j>! . (t) evaluated at t = -1, remains unchanged for 

1 > J 

all j. Each entry of the initial array ^ J is the value of the initial 
approximation to <J>! (t) evaluated at t v , . Specifically, according to the 

A 1 JL 

way algorithm I was constructed, this array has the elements 


*i,k,i = V t k-i' ri i>V 


’V 


(i = 1,2, 


,n; k = 1,2, 


,M+2) 


associated with the initial approximation 

$ 2 . q (t) “ (1 = 1,2, * • • »a) 


However, from the point of view of convergence, it may be advantageous to 

interrupt periodically the computation of the iterative procedure and restart 

it by supplying to the algorithm a new initiating array [4). , ]. The con- 

1, K, l 

struction of Steffensen's sequence provides a clue as to how this process of 
iteration modification can be applied effectively. An approach to the solu- 
tion of this problem might be attempted by first forming various sets of 
sequences corresponding to different groupings into disjoint subsets of all 
the entries of the array [i|^ j • In particular, for ease of implementation 

we consider the four sets of sequences consisting of 


(a) 


The corresponding elements 



(j = 1,2,, 


,n; k = 1,2, . . , ,M+2) 


(b) The total array |Vi,j,k 

(c) The corresponding rows 


If 


h, 


1, j * 


A 


i,M+2 




(d) The corresponding columns 


(i = 1,2, . . . ,n) 


(59) 


(60) 


i(*l,k,j**2,k,j» * • • »V,k,j) T |. (k = 1,2, , , . ,M+2) (61) 

3 = 0 

If certain assumptions, which may or may not be true, are made about each 
of these sequences, the initial array ^ jl can be modified by any one of 
four different ways. ■ *' * -* 
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Modification of individual entries .- If each sequence of equation (59) 
can be assumed to be iterates of a single-point scalar iteration function, 
each element of the initiating array [ij>. . ] can be modified by Aitken's 

formula (53) one at a time. Thus, 1 * * 2 

U. „ X . ♦. „ - ... 

new i,k, 3 


*i,k,3 - 2 *i,k,2 * 


i.k,'- 


For cases n > 2, however, numerical results show that such a procedure not 
only failed to hasten convergence but caused the algorithm to diverge. This 
is not surprising because too much information is lost due to our disregarding 

implicitly the fact that each element of an array of jkjJ j_ 0 is a func- 

tion of all the elements of the preceding array. The case with n = 1 is a 
special case of modification of columns to be discussed in a later paragraph. 


Modification of the total array .- In contrast to the method of modifica- 
tion of the last paragraph, the entire rectangular array can be con- 

sidered as an iterate of a single-point vector iteration function of n(M + 1) 
components. Then the initiating array [4'i ) j > k3 can be me^i-fied in toto by 

means of the matrix equation (57). This fully takes into account the depen- 
dence property mentioned in the last paragraph. Nevertheless, such a scheme 
could hardly be considered feasible from the standpoint of efficient computa- 
tion for the following reasons. Even if n and M are moderate in size, the 
matrix to be inverted is of a high order n(M + 1). Furthermore, the first 
modification cannot be effected until after n(M + 1) + 1 iterations by which 
time the algorithm has already converged or will have converged to a solution 
after a few more iterations. The case for n = 1 is more tenable but it sim- 
plifies to the method of modification of rows to be discussed presently. (The 
abbreviations MR and MC in the sequel denote the methods of modification of 
rows and modification of columns, respectively.) 

In view of what is said in the preceding paragraph, any acceleration of 
convergence procedure involving either modification of individual entries or 
the total array must be ruled out. On the other hand, the methods MR and MC 
intermediate in complexity between the two previous ones mentioned are found 
to be effective in the acceleration of convergence of algorithm I as well as 
algorithm II of the following section. 

Modification of rows (MR) .- Suppose each of the sequences (60) is an 
iterate of a single-point vector iteration function of M + 1 components. 

Then the initiating array [i/>. , ] can be modified one row at a time as 

follows: 1,5 


Let [<K ..] (j - 2,3, . . . ,M+3) be the M + 2 rectangular arrays 

1 > 9 J 

generated by the algorithm for an initiating array ^ j] . For each i, set 


y = ($ • i a • o • 

y j-i , 2,3 


»ip 


i,M+l , 


Cj = 1.2, 


,M+3) 
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fx*+ \ ^ 

Compute yo ' by means of the matrix equation (57) . The ith row of the 
new initiating array is then formed by setting 

- Xk" 1 ’ » * 1.2. • • • .M+l) 

After all the rows of [ij^ ^ x ] have been modified, the iterative process is 

restarted by returning to step 5 of algorithm I, The testing at step 8 should 
be bypassed until ^ 2 ] has been generated in the following iteration. 

Numerical results suggest that unless the sequence of rectangular arrays 
are already quadratically convergent within M + 2 iterations with no modifica- 
tion, MR modification, in many instances, substantially reduces the number of 
iterations required for convergence. It is rare when modification is needed 
more than once. 

Since the methods MR and total array modification coincide for n = 1, 
the latter method can be considered as a special case of the former in this 
instance. 


The feasibility of MR as a convergence acceleration procedure lies in 
the fact that it does not have the complicating features of total array modi- 
fication, and unlike entry-wise modification it evidently does make sufficient 
use of the history of previous iterations to be effective in the acceleration 
of convergence. 

The following method provides another and possibly more effective 
procedure in the acceleration of convergence for both algorithms I and II. 

Modification of columns (MC) .- Since MC modification is also inter- 
mediate in complexity between the first two methods discussed, it possesses 
the desirable features of MR modification. Moreover, when n is smaller 
than M, as often is the case in practice, it may be preferable to modify the 
initiating array -J by the former rather than the latter method. The 

validity of the statement will become evident in the ensuing discussion. 


The assumption made here is that each of the M + 1 sequences (61) are 
generated by a single-point vector iteration function of n components. Mod- 
ification in this case is carried out after every n + 1 iterations by making 
use of the n + 2 rectangular arrays [^j^k] (j = 1,2, . . . ,n+2) saved in 

storage. Each column (specifically the kth) of the initiating array is 
re-evaluated by setting 


(Vk.r^.kj 


■Vk.j)’ 


Cj = 1,2, 


,n+2) 


and determining y^ r+1 ^ by means of the matrix equation (57). The column in 
question is then replaced with new entries by setting 
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(i — 1,2, • , • , n) 




i,k,l 


.(*+ 1 ) 

i,o 


The iterative process, after every column of [>• v . ] has been thus modified, 

X , R. , 1 

as in the case of MR modification, is restarted at step 5 of algorithm I. 


The remarks made for MR modification about the testing at step 8 is 
also applicable here. 


Since M is usually set larger than n in practice, to achieve a 
desired accuracy, the MC method, besides saving storage allows modification 
to take place sooner and more frequently than modification by the MR method. 
When this occurs, numerical results show that convergence is accelerated more 
strongly than by MR modification. 


A comparison of the rates of convergence (in terms of number of itera- 
tions and machine time applied to the solution of example I) for iteration 
with MR and MC modification as well as no modification is given in table IV. 
It shows that both MR and MC methods substantially reduced the number of 
iterations and machine time required for convergence. It also shows that both 
the total number of iterations and the total machine time over all of the 
indicated intervals for MC modification is significantly less than that of 
MR. This should not, however, mislead us into precluding MR modification as 
a tool in the acceleration of convergence. Although numerical results show 
that whenever both methods work, MC modification works better; it is quite 
conceivable that acceleration of convergence in certain cases may fail for MC 
modification and work well for MR modification. For specific examples a 
proper combination of methods (with or without modification) may optimize 
the computation time required. 


TABLE IV.- A COMPARISON OF CONVERGENCE FOR ITERATION WITH AND WITHOUT 


MODIFICATION (ALGORITHM I) ■ 


j • , c -~'~Tr^*Interval 
Mo di f l c at l oh'-i-LL^^ 

3 < x < 4 

4 < x < 5 

5 < x < 6 

Total 

None 

38(5.16) 

44(6.00) 

53(7.32) 

135(18.48) 

MC 

19(3.42) 

18(3.30) 

14(2.52) 

51(9.24) 

MR 

25(4.50) 

15(3.06) 

15(3.12) 

55(10.68) 


Each entry of the table indicates the number of iterations required by 
algorithm I to solve example I for the indicated method of modification and 
interval of integration. The machine time in seconds on the IBM 7094 for the 
corresponding iteration is given in parentheses. The degree of approximation 
is M = 11 and the prescribed convergence error is e = 0.5xl0 -1 ^. 
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An Alternate Procedure 


Instead of using equations (42) and (43) , it is sometimes preferable from 
the point of view of faster convergence to construct an algorithm to solve a 
system of n first-order differential equations based on the equations 


■VP** 


■ 


- " 2 * J *2<“>h>l'h,r 




y (62) 


+ ” Tln+ J_ 1 , j + 1 , j + 1 * * * ‘ ^n-ij + l’^n.j 


)du 


with 


= r\ ± (i = 1,2, , . . ,n) 


(43) 


These equations differ from equations (42) and (43) in that the most 
recent information is used at each step. It will be shown in the next section 
how the specialization of (62) serves as the basis for a more efficient 
algorithm to solve an nth-order differential equation. 


INTEGRATION OF AN nth -ORDER DIFFERENTIAL EQUATION 


An nth-order differential equation can be expressed as a system of n 
first-order differential equations. Consequently, numerical methods for 
solving a system of first-order equations can be employed to solve an nth- 
order differential equation (see, e.g., example I). However, in so doing we 
may be performing many more operations which are otherwise unnecessary if the 
solution of nth-order differential equation were found without changing it to 
a system. This is true at least for the case of the Chebyshev series and will 
be explained in the ensuing discussion. 

Let the nth-order differential equation to be solved assume the form 

* (n) (t) = *£,*,♦! ... ,*("<>) • - lst£l (63) 

with the initial conditions 

<j> (l) (-l) = m (i = 0,1, . . . ,n-l) (64) 
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To provide the basis for the construction of the algorithm for the 
solution of an nth-order differential equation, define 


<h(t) = 

* 2 .(t) = *Cn-2)(t) 


<i> n CtO = *(t) 

Then equations (62) and (43) specialize at once to the equations 


(n-i) 
1 + 


(t) 


*' n ' 2) m 

j+1 



» 


*n 


du 




(65) 


4 


j+l 


(t) = 





with 


* 


(i) 

o 



(i — 0,1, . . . , n-1) 


( 66 ) 


The sequence of functions { ^ (t) } (as its counterpart 

l J ' j =o 



of 

j=0 


eqs. (42) and (43) ) is expected to converge uniformly to a function <f>(t) sat- 
isfying equations (63) and (64) on the closed interval [-1,1], Note that 
unlike the case where an nth-order differential equation is treated as a sys- 
tem of n first-order differential equations, the integration involved in 
each of equations (65) except the first is applied to the derivatives of the 
current iteration. Numerical results show that this gives algorithm II a 
clear advantage over algorithm I in terms of the number of iterations required 
to solve an nth-order differential equation for a given convergence error 
(see example II). Furthermore, if the Chebyshev series of 


^(t) = i|^t, <{>.,<}>:, 
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is available, the n integrals of (65) can be obtained readily by making use 
of the relation (B19) between the Chebyshev coefficients of the function and 
those of its derivatives. Hence, if an algorithm similar to algorithm I is 
constructed for the nth-order differential equation, the modified interpo- 
lating polynomial is used only to approximate the nth derivative during one 
iteration. This is in contrast to algorithm I which utilizes Qjq (t ) n times 
in the approximation of first derivatives. Since the computation of the modi- 
fied interpolating polynomial takes up the bulk of the time per iteration, the 
difference in time spent solving the same differential equation may be 
considerable. 


Algorithm II (for an nth-order equation) 


The construction of the algorithm for an approximate solution of an nth- 
order differential equation, as in the case of algorithm I, is based on the 


(i = 0,1, . . . ,n) 


j = 0 


assumption that the sequences of functions |$j^ (t)j 
can be represented accurately by polynomials of sufficiently high degree. Let 

M+n-i 


k=0 

denote, respectively. 


M+n- i 

Z) b k i)r k'« a " d 

k=o 


(67) 


.(i) 


.U) 


and *£J(t) 


(i— 0,1, . . . , n) 


( 68 ) 


For notational simplicity, equality between the respective expressions (67) 
and (68) is assumed to hold. The approximate solution of equation (63) with 
initial conditions (64) for a prescribed convergence error e is obtained as 
follows : 

1. Set 

k = 0, 1, . . . , M + 1 
points where ^ as an extremura -) 

2, Set 

b Q (i) = 2n ± (i = 0,1, . . . ,n-l) 

b^ =0 (i = 0,1, . . . ,n-l; k = 1,2, . . . , M+n-i) 

(This is equivalent to the initial approximations <}>^ 1J (t) = 

(i = 0,1, . . . , n-i).) 


. kTT 

h = cos mTT » 

(This computes the M + 2 
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3 . Compute 
M+n-i 


f J ct*) - E’ tfV*, 

r=o 


) (i = 0,1, . . . ,n-l; k = 0,1, 


. ,M+1) 


4. Compute 

♦£ fa • 

(k = 0,1, . . . ,M+1) 

5. Compute (using eq. (B27)) 

M+l , .. 

B cn) = 2 Y” ^\(t r n x (t k ) 

k M+l JL*/ J 




(k =0,1, . . . ,M) 


r=o 


(n) 


(We have here approximated <f>j + ^(t) by the Mth degree polynomial 

M , 


where according to theorem B6 
l 


B (n) 

B k * 


| J • • • ,f (n_i:) (t)]T k Ct)(l - t 2 )' 1 / 2 dt 


for a sufficiently large M.) 

6. Compute (by means of theorem B4 in appendix B and eq. (A31)) for 
i = n-l,n-2, . , . ,0 


g( i+ i) 

B (i) _ M+n-i-i 
M+n-i 2(M + n - i) 



R (i+1) R (i+l) 
B k-l - B k+l 


2 k 


(k = 1,2, . . .,M + n - i - 1) 
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(We have in this step obtained the n integrals 
M+n-i _ M 

' B k i)T k (t) ° f ^2 B k T k (t) (i = °’ 1, • * * >n ~^ 
k=o k=o 



7. If - b^ n ^ < e for all k, we are through. 

8. Otherwise, set 

b k X) = B k X) (i = °> 1 > * « * > n ; k = M, • • • , M+n-i) 

and return to step 3, 

Note: To take into consideration the case when 4 ^ (t) = (t) = 0 step 7 

should be skipped until the second iteration. 

Example II 

To provide a basis for comparison for both algorithms, we shall solve the 
differential equation (48) by algorithm II for the same prescribed convergence 
error and subintervals as in example I (i.e., e = O.SxlO -10 and xq = 0, Xj = 1, 
. . . , x 6 = 6). 

For each subinterval [x r _ 1 ,x r ] (r = 0,1, . . . ,6), let x = ct + d, 
c = y(x r+1 - x r ), d = j(x r + x r+1 ). Thus, 

F (l) (x) = F (l) (ct + d) = 2 1 <t> Cl - ) (t) (i = 0,1, 2, 3) 

Substitution into (48) then yields 

<r'Ct) = | (tret)] 2 - 2<j) (t) <f>" (t) - 

the form required by the algorithm. The solution is found as six initial- 
value problems corresponding to the six subintervals. The initial conditions 
used for each subinterval are determined by the solution at the end point of 
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the previous subinterval except the first where 


4> Ci) (-l) = \ F (i) (0) (i = 0,1, . . , ,n-l) 

2 1 

A comparison of computer results obtained by algorithm II with that of 
algorithm I shows agreement to 10 decimal places. Computer results also show 

that the coefficients B^^ (i - 0,1, . . . ,n) exhibit no change to 10 

decimal places for M > 11. Consequently, from the discussion in the 
preliminary analysis section the approximate Chebyshev series satisfies the 
differential equation with an error bound of order e = 0 . 5*10" 10 . 

A comparison of tables III and V reveals that the number of iterations 
required for the solution of the aforementioned differential equation is 
consistently less for algorithm II than for algorithm I. 


TABLE V.- THE INFLUENCE OF M ON THE CONVERGENCE OF ALGORITHM II 


^^Interval 

M 

0 < x < 1 

1 < X < 2 

2 < x < 3 

3 < x < 4 

4 < x < 5 

5 < x < 6 

11 

10 

17 ' 

24 

33 

39 

40 

12 

10 

17 

23 

31 

37 

37 

13 

10 

17 

23 

28 

34 

27 

14 

10 

17 

23 

26 

27 

26 

15 

10 

17 

23 

27 

28 

26 


Each entry of this table indicates the number of iterations algorithm II 
required to solve equation (48) having the indicated M and interval of inte- 
gration. The convergence criterion is e = 0.5xl0“ 10 . Since the approximate 
Chebyshev coefficients agree to 10 decimal places for all M > 11, the 
accuracy of the solution for each M is good to 10 decimal places. 


Remarks on Convergence 
Since a polynomial is used to approximate 

<t>|^(t) = <J^t, $.,((>!, . . . 

the convergence of algorithm II, as in the case of algorithm I, may be influ- 
enced by the degree of the approximating polynomial M. Examination of 
table V shows that while the number of iterations needed for convergence for 
the first three intervals remains essentially unchanged, the number of iter- 
ations in the remaining three intervals required for convergence is essen- 
tially a decreasing function of M. Thus by a judicious choice of M the 
computing time spent for a particular problem may be minimized. However, 
when such a choice is unavailable, it is probably better to pick a higher M 
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than is thought necessary so that convergence can be achieved within a toler- 
able number of iterations. For example, in the case of example II, it would 
have been better to pick M = 15 instead of 11. 

It will be shown next that the convergence of algorithm II for a fixed 
M can be accelerated by methods already proposed for a system of n first- 
order differential equations. 


Acceleration of Convergence 


Note first that the sequences 4>j n ^ 


(t k )| (k = 0,1, . . . ,M + 1) 


can be considered as a single sequence of one-dimensional array 




» 4 > 


j ,M+2^ (j - 1j2, . . . ) 


(69) 


where 


ip . , = 4> (t, ) 

,k v k- r 


(70) 


Thus if one thinks of (69) as a rectangular array of one row and M + 2 
columns, the methods of MR and MC modification proposed for algorithm I can 
also be used in the acceleration of convergence for the nth-order differen- 
tial equation. The initiating array 1 ,\p 2 j, . . . ^ is , 

re-evaluated for every M + 2 or every two iterations depending upon whether 
MR or MC modification is utilized. Since MR modification is applied to a 
single row and MC is actually a modification of individual entries here, the 
number of operations involved compare with that of algorithm I (applied to the 
solution of the same nth-order differential equation) should be considerably 
less. 


Table VI shows that the number of iterations as well as the machine time 
required in the solution of example II for both MR and MC methods is consid- 
erably less than that of iteration without modification. It also shows that 
MC modification has a definite edge over MR modification both in the number 
of iterations and in machine time required for convergence. However, the 
comments applied to algorithm I about not precluding MR modification as a 
tool in the acceleration of convergence should also be heeded in the present 
case. 


Table VII is based on the data of tables IV and table VI and illustrates 
the advantage of algorithm II over algorithm I in the amount of machine time 
spent on solving the same differential equation. It shows that in every case, 
whether modification is involved or not, the machine time required by 
algorithm I is at least 1-1/2 times that of algorithm II. In the case of MC 
modification the ratio of the machine time of algorithm I to that of 
algorithm II is as large as 2.86. 
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TABLE VI.- A COMPARISON OF THE CONVERGENCE FOR ITERATION WITH 


AND WITHOUT MODIFICATION (ALGORITHM II) 


.. j .T^7i^~-~-Jnterval 
Mod i f i c at 

3 < x < 4 

4 < x < 5 

5 < x < 6 

Total 

None 

33(3.25) 

39(3.86) 

40(3.98) 

112(11.09) 

MC 

16(1.80) 

14(1.57) 

8(0.88) 

38(4.25) 

MR 

18(2.21) 

15(1.92) 

15(1.93) 

48(6.05) 


Each entry of this table indicates the number of iterations required when 
using algorithm II to solve example II for the indicated method of modifica- 
tion and interval of integration. The machine time in seconds on the IBM 7094 
for the corresponding iteration is given in parentheses. The degree of 
approximation is M = 11 and the convergence criterion is e = 0.5xl0~ 10 . 


TABLE VII.- RATIOS OF MACHINE TIME OF ALGORITHM I TO ALGORITHM II 



Each entry of this table gives the ratio of the machine time of 
algorithm I to that of algorithm II based on tables IV and VI. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, May 15, 1969 
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APPENDIX A 


CHEBYSHEV POLYNOMIALS 


This appendix is included to facilitate the discussion of the 
approximation of functions in terms of Chebyshev polynomials. It also pro- 
vides the minimum of the necessary working tools for most numerical work 
involving these polynomials. Most of the properties can be found in the works 
of Rivlin (ref. 10) and Lanczos (ref. 11) . Clenshaw's summation formula for 
Chebyshev series is generalized here for the summation of sequences obeying an 
nth-order linear recurrence relation (see theorem Al) . 


Some Properties of Chebyshev Polynomials 

Definition . - The polynomial of degree k defined by 

T^(t) = cos (k arc cos t) , -1 < t < 1 (Al) 

is called the Chebyshev polynomial of the first kind of order k or simply a 
Chebyshev polynomial. 

As an immediate consequence of this definition we have. 

Property Al . If k * 2, then 

T k (t) = 2tT k _ x (t) - T k _ 2 (t) (A2) 


with 


T 0 (t) = 1, T r (t) = t 


From the recurrence relation (A2) we can generate Chebyshev polynomials 
up to any order n by setting k = 2, 3, . . . , n. 

By the change of variable t = cos 0 one can demonstrate that the 

sequence of Chebyshev polynomials I Tv. (t) } is orthogonal on [-1,1] with 

l K I k=0 

respect to the weight function w(t) = (1-t 2 ) -1 /2 , that is, 

Property A2 . For any two Chebyshev polynomials T .(t) and T k (t')» the 
following conditions hold. 


P T m (t)T k (t)(l - t 2 )' 1 / 2 dt = 


if m ^ k 

if m = k ^ 0 

if m = k = 0 


(A3) 
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Consequently, any results deduced for general classes of orthogonal 
polynomials also hold for Chebyshev polynomials. 

Property A3 (Orthogonality with Respect to Summation) , 


Let 



cos 


3JL. 

N + 1 


(j = 0,1, . . , ,N+1) 


(A4) 


If T m (t) and T k (t) are any pair of Chebyshev polynomials of orders zero 
through N + 1, then 

m = k = 0 or N+l 
m = k ^ 0 or N + l (A5) 

m ^ k 

Furthermore, if T r (t) is a Chebyshev polynomial of any order 


N+l 


VWV - 


J rrf) 


N + 
(N + 
0 


1 

l)/2 


if 

if 

if 


E ' 


Tr(tl) 


j=0 


N+l if r = 2s (N + 1) , 
0 if r / 2s (N + 1) , 


s - 0,1,2, 
s — 0,1,2, 


(A6) 


Property A4 . For all k > 1, T^Ct) has zeros at 


t , cos C2J * 1)» 

t cos 2k 

(j = 0,1, , . , ,k-l) 

(A7) 

and extrema (-l)* 1 at 


(j = 0,1, . . . ,k) 

(AS) 


Property A3 together with equation (A8) is instrumental to the discussion 
of the "Modified Interpolation" of appendix B. The following property is use- 
ful in the integration of a function represented by a series of Chebyshev 
polynomials: 

Property A5 . 


jT k Ct)dt 


Tj(t) 

T 2 (t)/4 


1/2 


y. k + 1 



k = 0 
k = 1 


k > 2 


(A9) 
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Summation of Series of Chebyshev Polynomials 


We consider next the summation of a finite series of Chebyshev 
polynomials of the form 

N 


- S VkCt) 

k=o 


(A10) 


where Bj, are constants. The evaluation of S^(t) can be efficiently 
performed by particularization of the following theorem (see also ref. 12): 
Theorem Al. If 


N 

= 2 a k u l 


k=o 

where u^ obeys the nth-order recurrence relation 

n 

u k » k - n 

3 = 1 


(All) 


(A12) 


then S can be evaluated by the formula 


s = c o u o + yi c kf u k - y^ r j,k u k_i\ 
k=i V j=i / 


(A13) 


where c Q ,c 1; 


,c n _ 1 are determined from the recurrence formula 


c N+l “ c N +2 c N+n = 0 


n 

c k = a k + ^ > r j,k+j c k+j 
3=1 


(k = N,N-1, . . . ,0) 


(A14) 


The theorem can be established as follows: Replace a^ in equation (All) 

by equation (A14) so that : 


N / 3* 

■EK-E-j, 


k+j c k+jl u k 


k=0 


3=1 


or 
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s = 


N n 


E ** -E 

k=o j=l 


XI r j ,k+j c k+j u k 

k=o 


Setting now k + j = m in the double sura and noting that = 0 for 
k=N+l, N+2, . . . , N + n, we have 


S = 


N n 

E <** - E 


k=0 j=i 


N 

X! r j ,m c m u m- j 

m=j 


or by inverting the order of summation for the double sum, one obtains 
N iH m N n 

s - E c « E E r j ,» c mVj * E E r s 

k=o m=i j=i m=n j = i 

The collecting of terms in c m then yields 


S = c 


n-i / v N / n \ 

> u o + E K E r j.«Vjk +EK-E r i,. u .-j) 

m=i \ j = i / m=n \ j = i / 


m 


(A15) 


But by equation (A12) 

n 

"m -2j r j,m u ra-j = 0 f°r m > n 

j=i 

This proves that equation (A13) holds. 

Useful special cases : 

1. For n = 1, we have 


u k = r k u k-i < k * 1) 


CA16) 


so S of (A13) becomes 

S = c 0 u Q CA17) 

with 

C N+1 = 0 > c k = a k + r k+l c k+l (k = N,N-1, ... ,0) (A18) 


39 



2. For n = 2, we have 


U k = r i,k U k-l + r 2,k U k-2 

so S of equation (A13) becomes 

S = c Q u 0 + Cl ( Ul - r 1}1 u 0 ) 

with 


C N+ 1 c N+2 " 0 


c k = a k + r i,k+i c k+i + r 2,k+2 c k+2 
(k = N,N-1, ... ,0) 




The following properties are applications of the special case of 
to sums of Chebyshev polynomials: 


Property A6 . If 


N 

s(t) =2 B k T k Ct) 


then also 


k=o 


S(t) = ~ + c : t - c 2 


where Cj and c 2 are generated by the recurrence relations 


'N+i = c N+2 " 0 


c k - B k + 2tc k+1 - c k+ ^ (k - N,N-1 , . . . ,1) 

(This formula follows from eq. (A21) and the fact that 

T k (t) = 2tT Jc _ 1 (t) - T k _ 2 (t) (k > 2)) 

Property A7 . If 

S(t) - 2’ B k T 2k'« 


k=o 


(A19) 


(A20) 


(A21) 
= 2 


(A22) 


(A23) 


(A24) 


(A25) 
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then also 


B n 

SCt) = -f + c x (2t 2 - 1) - c 2 

where c^ and c 2 are generated by the recurrence relations 

C N+ 1 = c N+2 = 0 | 

c k = ® k + ^ (2t 2 - 1) c k+ ^ - c k+2 Ck = N,N-1, . . . , 1) | 

(The formula follows from eq. (A21) and the fact that 

T 2k (t) = 2 (2t 2 - l)T 2k . 2 - T 2k . 2 (t) - T 2k _ t (t) (k > 2)) 

Property A8 . If 

SCt) = ^B k T 2k+1 (t) 
k=l 

then also 

S (t) = (2t 2 - l)( Cl - c 2 ) 

where c x and c 2 are generated by the recurrence relations 

> 

C N+1 = c N+2 = 0 

► 

C k = Bk + 2(2t 2 - l)c k+1 - c k+2 (k = N,N-1, ... ,1) 

> 

(The formula follows from equation (A21) and the fact that 

T 2k*l<« ' 2 < 2t2 - DTzk-lC*) - T 2k-3tt)) 

Formulas for evaluating S^ft) of equation (A10) at special values 
are given as follows: 

Property A9 , If 

N , 

%(t) = B k T k (t) , -1 < t < 1 

k=0 


(A26) 


(A27) 


(A28) 


(A29) 


(A30) 


Of t 


(A10) 
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then 


and 


k=o 


%C0) = 2' c-^S 


%CD = ]>j 


(A31) 


(A3 2) 


(A33) 
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APPENDIX B 


THE CHEBYSHEV SERIES AND APPROXIMATION BY MODIFIED INTERPOLATION 


The Best Approximating Polynomial and the Chebyshev Series 

Let f(x) be a continuous function defined on a closed interval [a,b] and 
let e be a prescribed positive number. The existence of a polynomial P(x) 
for which 

max |P(x) - fCx)| < e (Bl) 

a<x<b 

is given by a well-known theorem of Weierstrass. We introduce here the con- 
cept of the best approximating polynomial to facilitate later discussion. 

(For a detailed discussion of the best approximating polynomial and its prop- 
erties, see ref. 13.) 

Definition .- Let D^ denote the set of polynomials of degree s N. A 
polynomial P*(x)eD N having the maximum residual: 

max | f (x) - P*(x) j = min / max^ jf(x) - P(x)A (B2) 

a<x<b 1 1 PeD n \a<x<b 1 ') 

is called a best a-pprox-mating polynomial of degree N of f(x) in the 
Chebyshev sense. 

Definition .- The quantity 

EN(f) = ?IB n (a?gb |f(X) ' P(x:)| ) CB3) 

is called the smallest deviation of the polynomials of D n from f(x) or the 
minimax. 

The best approximating polynomial always exists and is unique. It is 
completely characterized by the theorem of P. L. Chebyshev. 

Theorem Bl . Let f(x) be a function continuous on [a,b] . Then any 
polynomial P(x)eD N is the best approximating polynomial if and only if 
there exist N + 2 points 

a < x 1 < x 2 < . . . < x N+2 < b (B4) 

for which 

|f(xi) - P(xi)| = max |f(x) - P(x)| = E (i = 1,2, . . . ,N+2) (B5) 

a<x<b 

with the function f(x) - P(x) alternating in sign at consecutive values of 


A useful consequence of the above theorem is: 


43 



Corollary . - If f(x) is continuous on [a,b] and if for any Q(x) in D N 

the function f(x) - Q(x) alternates in sign on a set of N + 2 distinct 
points 


a < Xj < x. < 


< 


X N+2 


< b 


CB6) 


with 

|f Cx ± ) - QCx ± ) | = M (B7) 

then 


M < E N (f) 


(B8) 


Let ij)(t) be continuous on the closed interval [-1,1]. We shall be 
interested in the expansion 

00 

<Kt) = ^ A k T k (t) , -1 < t < 1 (B9) 

k=o 


where T k (t) are Chebyshev polynomials (see appendix A) defined by 

T k (t) = eos(k arc cos t) (BIO) 

Definition . - In the particular case where 


A k = | 4>(t)T k (t) (1 - 1 2 )' 1/2 dt (Bll) 

the series (B9) is known as the Chebyshev series for $ (t) , and the coeffi- 
cients A k are called Chebyshev coefficients . (See refs. 10 and 13 for a 

detailed treatment of expansion in Chebyshev series, and of certain relations 
between the best approximating polynomial and the Chebyshev series.) 

The following exhibits a large class of functions which have uniformly 
convergent Chebyshev series expansions: 

Theorem B2 .- If (t) satisfies the Holder condition, that is, if there 
exists constants M and a such that, for all tj and %2 in [-1,1] 

UCtD - <f»(t 2 )| < M|t! - t 2 | a (a > 0) (B12) 

then 4> (t) can be expanded in a uniformly convergent Chebyshev series. (For 
efficient computation, however, the given function should have stronger prop- 
erties such as differentiability. Consequently, prior to expansion in 
Chebyshev series, it may be necessary to provide such properties by suitable 
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transformations or subdivision of the interval of definition.) 

Let the partial sum formed by the first N + 1 terms of the Chebyshev 
series of a function 4>(t) be denoted by 

N 

= E' C B13 > 

k=o 

and the maximum error of S^(t) by 

a N C<j>) = max |<Kt) - S N (t) | (B14) 

-l<t<l 

we note here that, since S N (t) is a linear combination of polynomials at 
most degree N, the partial sum S N (t)eD^. 

The following theorem gives an important inequality between the minimax 
error E N (4>) and the error a N (<j>): 

Theorem B3 (A. Lesbesque) .- If <j>(t) is continuous on [-1,1], then 

a N (<)0 < (3 + log N)E n O) (B15) 

The above inequality means that for practical purposes the truncated 
Chebyshev series is just as good as the best approximating polynomial. 

Some useful inequalities for a function <J)(t) which is expandable in a 
uniformly convergent Chebyshev series are: 


e n O) < 






(B16) 


max | P* ft) 
•l<t<l 


s N (t)| - 2 k 2-f l I A : 


/ oo w 

1/2 E s s E iv 

k=N+l k=N+l 


(B17) 


(B18) 


A useful relation between the Chebyshev 
the Chebyshev coefficients of its derivative 


coefficients of a function and 
is given by: 
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Theorem B4 .- Let <(>(t) be defined on [-1,1]. If <j>’(t) is integrable, 

then 

A Cl) „ A (l) = 2kA iC0) (k =1,2, . . .) (B19) 

where A^ 0 -* and are Chebyshev coefficients of <f>(t) and <f>'(t), 

respectively. 

(The validity of eq. (B19) can be demonstrated by making the change of 
variable t = cos 0 in eq. (Bll) and then integrating by parts.) 


Approximation by Modified Interpolation 

Let f(x) be continuous on [a,b] and let denote the set of all 

polynomials of at most, degree N. Recall that the best approximating poly- 
nomial of f(x) as defined in the preceding section is the polynomial 
P*(x)eDn for which 

max |f(x) - P* (x) I = min ( max, |f(x) - P(x)|) 
a<x<b 1 1 PeD N \a<x<b 1 / 

It has also been stated that the necessary and sufficient condition for a 
polynomial P(x)eD^ to be a best approximating polynomial is that there 

exist in [a,b] at least N + 2 points 

x i < x 2 < . . . < x n+2 

for which 

|f(x.) - P(Xi)| = max | f (x) - P(x,)| = E (i = 1,2, . . . ,N+2) 

1 x a<x<b 1 

with f(x) - P(x) alternating in sign at consecutive values of x^. In view 
of this and the fact that T^ +1 (t) assumes its extrema (-1)^ at the points 

t j = cos WTT » j = 0,1 > • ' * ,N+1 

it is easy to see that the following holds: 

Lemma.- If 


N+l 

W* 5 “D'Vkt 1 ) 

k=0 


(-1 < t < 1) 


(B20) 
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then the polynomial 


N 


V 1 ’ 


r 

k=o 


Bk^VC^) 


(B21) 


is the best approximating polynomial of at most, degree 
the function B^ +1 t n+i W assumes its extrema (-1)‘ , 3 N+1 


N to G N+1 (t), and 
at the points 


t. 

J 


cos 


iz 

N + 1 


Cj = 0,1,2, . . . ,N+1) 


(B22) 


The interpolating polynomial P^ +1 (t).- Attention is now turned to the 

interpolating polynomial used by Clenshaw in reference 3. A modified version 
of it is a basic tool of the algorithms of this report. His polynomial 
assumes the form 


N+l 

w*> 

k=o 




(B23) 


and interpolates a given function <ji.(t) at N + 2 points of t, given by 
equation (B22) . Thus the coefficients B k of P^j+i Ct) can be determined by 

solving directly the linear system 
N+l 

y) n B k T k Ct.) = Ktj) Cj = 0,1,2, . . . ,N+1) (B24) 

k=o 


of N + 2 equations in N + 2 unknowns . However, these equations can best 
be solved by taking advantage of the orthogonality property of Chebyshev 
polynomials with respect to summation. First multiply each one of the systems 
of equations (B24) by T m (tj) and then divide the first and last equations by 

2. Hence, upon adding the resulting equations, we have 


N+l N+l N+l 


j=o 


j=o k=o 


Changing the order of summation yields 

N+l N+l N+l 

2 J *C t j) T mCtp =£" B k 
j=0 k=0 j*0 
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and thus, by equation (A5), 

N+l 


B,. = 


k N+l 




(k = 0,1, . 


j = 0 


>N+1) 


(B25) 


Since 


W - cos r?r = W 


(B26) 


the coefficient Bj, can be written also as 

N+l 

TT ]C'' Kt PVV 


B, — , 
k N 


(k = 0,1, . . . ,N+1) 


1=0 


(B27) 


We note here that equation (B27) is more desirable than (B25) from the point 
of view of computation. Since B^ is a finite series of Chebyshev poly- 
nomials evaluated at tj., it can be readily calculated by means of Property A6. 

Some important properties of the approximating polynomial P n +,(t) are 
given in the following theorem. 

Theorem B5 .- If a function <f>(t) is continuous on [-1,1], and P N+1 (t) 
is the approximating polynomial defined by equations (B23) and (B27) , 
then for each k „ 


B k ~ A k + 


S ( A 2r(N+i)-k + A 2r(N+l)+k) 


(B28) 


r=i 


where A^ are the Chebyshev coefficients of 4>(t), and 


2 l%+l1 - E NC^ 


(B29) 


Equation (B28) can be demonstrated to hold by means of equations (A5) and 
(A6) . To prove (B29) note that by the lemma at the beginning of the section 


Q n (c) 


N , 

= 2 VkCt) 


k=o 


is the best approximating polynomial of maximum degree N to P^ +1 (t). The 
residual function P N+1 (t) - Q N (t ) = (l/2)B N+1 T N+l (t) assumes its extrema 
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(-l)-j (B n+ 1 /2) at the points t^ = cos[jTr/(N + 1)], (j = 0,1, . , . ,N+1). Now 

since 4>(t) coincides with P^ +1 (t) at the same N + 2 points, the residual 

function cf>(t) - Q^Ct) also assumes the values (-1)-* ( B N+1 /2) at these points. 

Recall that if f(x) is continuous on [a,b] and if for any Q^(x) i- n the 

residual function f (x) - Q^(x) alternates in sign on a set of distinct points 

a < x 1 < x 2 < . . . < x n+2 < b with ] f Cx ± ) - Q N Cx i ) | = M, then M < E^(f) . 

Hence it follows that (1/2}|B^ +1 | < E^(<|>). This proves equation (B29) and 
also the theorem. 

The modified interpolating polynomial Q^(t).- The polynomial P^ +i (t) is 
an interpolating polynomial having the N + 2 extrema of T^ +1 (t) as the 
points of interpolation. Let P^(t) denote the interpolating polynomial hav- 
ing N + 1 extrema of T^> (t) as the points of interpolation. We discuss next 

the modified interpolating polynomial Qjq formed by truncating the last 
term of P^ +1 (t); that is, " 

N 

(E30) 

k=o 

and also point out why Q^(t) is preferred over P^ft) as an approximating 
polynomial. 

Consider first the maximum deviation of P^ +1 (t) from S^ +1 (t), the first 
N + 2 terms of the Chebyshev series. Writing 


N+l 

Vi (t > - Vi") CBfc ■ Ak)TkW * ' Vi )T »*i t,) 

k=o 

and taking absolute values one obtains 

l P N+l^ " " Ak ^ + ^ B N+1 “ ^Sl+J 

k=0 


(B31) 


But we have by equation (B28) the inequalities 
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?I B 0 

- \1 

< 


1 A 2N+2 1 


+ 

1 A 4 N+ 4 1 

+ . . : 



I B 1 

- aJ 

< 

I A 2N+iI 

+ 1 A 2.N+3 1 

+ 1 A 4N+ 3 1 

+ 

1 A 4N+5 1 

+ * • . 




• 








\ 

(B32) 

I %-1 - 

a n-J 

< 

I A N+ 3 ! 

+ l A 3N+J 

+ t A 3N+5 1 

+ 

1 A 5N+3 1 

+ . . . 



l B N 

- a n | 

< 

l A N+2 

+ 1 A 3N+2 1 

1 A 3N+4 I 

+ 

1 A 5N+4 1 

+ - . * 



1 B 

2 B N+l ” 

a n+J 

< 


1 A 3N+3 1 


+ 

1 A 5n+5 1 

+ . . . 




Add the N + 2 equations and sum the right member according to ascending 
indexes to yield 

|Bfc - A kl + - A N.J ^k 1 (B53) 

k=0 k=N+2 

It follows from inequalities (B31) and B(33) that 

oo 

IW^ - W«l <- 53 |Akl (B34) 

k=N+2 

The application of the triangle inequality 

l<Ht) - P N+ 1 Ct)| < iKt) - s N+ 1 (t)| + |S N+I (t) - P N+ 1 (t)| 


together with inequalities (B17) and (B33) yields the error bound for P N+1 (t) 


max 

-l<t<l 


k(t) - p N+1 ct)| <2 V 

k=N+2 


|A k l 


(B35) 


We consider now the modified polynomial Q N (t) formed by the first N + 1 
terms of p jq +1 (t)- 1 


Since 


iQ N ct) - 8,(01*53 

k=0 


k - \i 
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it follows easily from (B32) that 


w 

max |<Kt) - Q N (t) | < a N (<f>) + max |Q N (t) - S N (t)| < |A n+1 | + 2 V' | A k | 


k=N+2 


On the other hand, we have by replacing N + 1 of inequality (B35) by N 


max 

-l<t<l 


OO CO 

IfKt) - P N Ct)| < ^2 |A k l < 2 |A N+1 | + Yj 


k=N+l 


+ ' l«kl 
k=N+ 2 


(B36) 


(B37) 


Hence it can be seen from inequalities (B36) and (B37) that the maximum error 
for P N Ct) for sufficiently large N can be two times larger than that of 

Qjyj(t). Moreover, since 


<Kt) - Q N Ct)j < |Kt) - P N+i Ct)| h- |P N+i (t) - Q N (t) | 


we have from the fact that l T N+1 (b)| < 1 and the inequality (B35) that 

oo 

;ax Uct) - Q N (t)| < \ |B N+1 | + 2 V |A k | (B38) 


k=N+2 


But 


2 J b n+i I S E N Cd>) < _ jfCt) - Q N (t) 


Consequently, 


oo 

\ l B N+ il S _ 1 5g 1 \Ut) - Q N Ct)j < \ |B n+1 I + 2 l A k l • t B39 ) 


k=N+2 


Note that since 


♦<v -r 


" 2 B N+l T N+l^j^ 


k=o 


by (24), the residual function <J)(t) - Q^(t) also assumes the values 
(-1)- 1 (B n+i / 2) at N + 2 points given by (B22) . 
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Thus if 


2 L i*k 

k-N+2 


is small relative to (1/2) |Bj, I (and this is often the case in practice), we 
have 1 

-l?t<l " Q N (t)l ~ \ t B N+ 1 i 

and Q n C t) closely approximates the best approximating polynomial of <Kt). 
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APPENDIX C 


FORTRAN IV SUBROUTINES FOR ALGORITHMS I AND II 


Subroutine AL ALG1 


Identification 


AL ALG1, Chebyshev Series Integration of a System of n First-Order Nonlinear 
Differential Equations 
FORTRAN IV, Double-Precision Subroutine 

Purpose 

This subroutine is used to generate an approximate Chebyshev series solution 
for a system of n first-order nonlinear differential equations with n ini- 
tial conditions. The differential equations are of the form 

d<J> • 

= • • • j = n) 

with the initial conditions 

= ri£ (i- = 1 j2, . . . ,n) 

The approximate Chebyshev series solution and derivatives are provided in the 
form of the finite series 

M+1-P 

^ p) (t) « B^T k Ct) (p = 0,1; i = 1,2, . . . ,n) 

k=o 

where T k (t) are Chebyshev polynomials defined by 

T^(t) = cos (k cos -1 t) , -1 < t < 1 

The accuracy of the approximate solution depends on the convergence error e 
and the degree of polynomial approximation M prescribed by the user. When 
both e -* 0 and M -> °°, the approximate Chebyshev series solution approaches 
that of the infinite Chebyshev series expansion (see main body of the report 
for choice of e and M and the estimation of errors). 
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Usage 


The routine is entered via the statement 

CALL ALG1 (N, KIN, M + 1, M + 2, ETA, EPSN, KIT, TR, PHI, XR, XS, B, IC, 
NER, DERIV) 


where 

N(=n) is the number of first-order differential equations. 

KIN is an integer code used to indicate the method of computation desired: 
KIN = 0, for straight iteration, (see algorithm I, page 13). 

KIN = 1, for iteration with modification of columns (see page 26). 

KIN = 2, for iteration with modification of rows (see page 25). 

M is the degree of polynomial approximation used to represent the deriva- 

tives <f^(t) = ^ i (t,c() 1 ,<j) 2 , . . . ,<f> n ), 

ETA is a double-precision array of 3n locations. The first n locations 
are used to store the n initial conditions (i = 1,2, . . . ,n) 

in the order of ascending i. The remaining locations are used 
internally by the subroutine. 

EPSN is the convergence error e prescribed by the user and is a double- 
precision variable. 

KIT is the maximum number of iterations allowed by the user. 

TR is a double-precision array of M + 2 locations reserved as working 

storage for the subroutine. 

PHI is a double-precision array reserved as working storage for the 
subroutine. The number of locations allocated are 


2n(M + 2) 


for 

KIN = 0 

n(n + 3) (M + 

2) 

for 

KIN - 1 

n(M + 4) (M + 

2) 

for 

KIN = 2 


XR is a double-precision array of n(4n + 2) locations reserved exclu- 
sively as working storage for the case with KIN =1. XR is a dummy 
double-precision variable for the cases with KIN = 0 and 2. 

XS is a double-precision array of (M + 1) (4M + 6) locations reserved 

exclusively as working storage for the case with KIN = 2. XS is a 
dummy double-precision variable for the cases with KIN = 0 and 1. 
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B is a double-precision three-dimensional array of dimension 

(M + 2) x N x 3. The first two rectangular arrays of dimension 

(M + 2) x N are used to store the approximate Chebyshev coefficients 

of <|>!(t) with 

b 5 0 . stored in B(K + 1, I, 1) 

K, 1 

for i = I = 1, 2, . . . , n; k = K = 0, 1, . . . , M + 1 

b/ 1 } stored in B(K + 1, I, 2) 

K, 1 

for i = I = 1, 2, . . . , n; k = K = 0, 1, . . . , M 

The remaining locations are used internally by the subroutine. 

IC is the number of iterations executed by the subroutine to achieve 
convergence. 

NER is the error code. NER = 0 is a normal return; NER = 1 indicates that 
the number of iterations had exceeded KIT. 

DERIV is the name of a user supplied subroutine for the computation of the n 
first derivatives <Jk (t) (see Derivative Subroutine below). The name 
DERIV (or whatever name the user chooses) must appear in an EXTERNAL 
specification statement in the calling program. 

Derivative Subroutine 


A subroutine for the computation of the first derivatives cf>!(t) 

(i = 1, 2, . . . , n) must be supplied by the user and must be of the 
following format: 

SUBROUTINE DERIV (N, T, PHI, PSI) 

DOUBLE PRECISION T, PHI, PSI 
DIMENSION PHI (N) , PSI(N) 

PSI(l) = . . . 

PSI (2) = . . . 


PSI (N) = . . . 

RETURN 

END 

The symbols of the DERIV subroutine are defined as follows: 

N(=n) is the number of first-order differential equations. 

T is the independent variable. 

PHI is an array of n locations used to store the values of (t) 
(i - 1, 2, . . . , n) in the order of ascending i. 
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PSI is an array of n locations used to store the n first derivatives 

<Ct) = if'iCt, $ 1 , <f> 2 > • • • > ♦n) Ci = 1» 2, . . . , n) in the order 
of ascending i. 

Other Subroutines Required 

The subroutines 

1. AL CHBY 

2. AL DPIN 

Evaluation of Solution and Derivatives 


The finite series 


M+l-p 

00 Bk!i T k Ct) ( P = °> 1 > 1 = 1 > 2 > • • • > n -> 

k=0 

may be evaluated as a function of the independent variable t by means of 
subroutine AL CHBY (see page 59), which in this case may be accessed via the 
statement 


CALL CHBY (B(l, I , P + 1) , M + 2 - P , T, SUM) 

where 


I, M, P 
B(l, I »P + 1) 
T 

SUM 


are integers set equal to i, M and p, respectively. 

is the location of B.^. 

1,0 

is the independent variable t. 
is the evaluated result. 


Subroutine AL ALG2 


Identification 


AL ALG2, Chebyshev Series Integration of an nth-order Nonlinear Differential 
Equation 

FORTRAN IV, Double-Precision Subroutine 
Purpose 

This subroutine is used to generate an approximate Chebyshev series solution 
for an nth-order nonlinear differential equation with n initial conditions. 
The nth-order differential equation is of the form 

♦ ^ (t) = *( • • - ,<t’ ( - n - i:) ) , -1 < t < 1 
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with the initial conditions 

<|> (l) (-l) = ni Ci = 0,1, . . , ,n-i) 

The approximate Chebyshev solution and its derivatives are provided in the 
form of the finite series 

M+n- i 

4 (l) (t) « 2 B k l)T k (t) Ci = °’ 1 ’ • ' • ’ n) 

k=o 


where T^(t) are Chebyshev polynomials defined by 

T^(t) = cos (k cos -1 t) , -1 < t < 1 

The accuracy of the approximate solution depends on the convergence error e 
and the degree of polynomial approximation M prescribed by the user. When 
both e -> 0 and M -*■ °°, the approximate Chebyshev series approaches that of the 
infinite Chebyshev series expansion. (See main body of the report for choice 
of e and M and the estimation of errors.) 

Usage 

The routine is entered via the statement 

CALL ALG2 (N, KIN, M + 1, M + 2, ETA, EPSN, KIT, TR, PHI, XR, B, IC, NER, 
NDER) 


where 

N(=n) is the order of the differential equation. 

KIN is an integer code used to indicate the method of computation desired: 
KIN = 0, for straight iteration (see algorithm II, page 30). 

KIN = 1, for iteration with modification of individual entries (see 
page 34). 

KIN = 2, for iteration with modification of rows (see page 34). 

M is the degree of polynomial approximation used to represent the nth 

derivative <j>^(t) = 4'(t,<|>><f>\> • . . ,4* . 

ETA is a double-precision array of 2n locations. The first n cells are 
used to store the n initial conditions rii (i = 0,1, . . . ,n-l) 
in the order of ascending i. The remaining cells are used inter- 
nally by the subroutine. 

is the convergence error e prescribed by the user and is a double- 
precision variable. 
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TR is a double-precision array of M + 2 cells reserved as working storage 
for the subroutine. 

PHI is a double-precision array reserved as working storage for the sub- 
routine. The number of locations allocated are 


(n + 

1) (M + 

2) 

for 

KIN = 0 

(n + 

3) (M + 

2) 

for 

KIN = 1 

(M + 

2) (M + 

n + 3) 

for 

KIN = 2 


XR is a double-precision array of (M + 1) (4M + 6) locations reserved 

exclusively as working storage for the ease with KIN =2. XR is a 
dummy double-precision variable for the cases with KIN = 0 and 1. 

B is a double-precision two-dimensional array of dimension 
(M + n + 1) x (n + 2) with 

B^ stored in B(K + 1, I + 1) 

for i = I = 0, 1, . . ... , n; k = K = 0, 1, . . . , M + n - i. 

The remaining cells are used internally by the subroutine. 

IC is the number of iterations executed by the subroutine to achieve 
convergence . 

NER is the error code. NER = 0 is a normal return; NER = 1 indicates that 
the number of iterations exceeded KIT. 

NDER is the name of a user supplied subroutine for the computation of the 

nth derivative <j> ^ (t) (see Derivative Subroutine below). The name 
NDER (or whatever name the user chooses) must appear in an EXTERNAL 
specification statement in the calling program. 

Derivative Subroutine 


A subroutine must be supplied by the user to compute the nth-derivative 
<|>^(t) and must be of the following format: 

SUBROUTINE NDER (N, T, PHI, PSI) 

DOUBLE PRECISION T, PHI, PSI 
DIMENSION PHI (N) 

PSI = . . . 

RETURN 

END 


The symbols of the NDER subroutine are defined as follows: 

N(=n) is the order of the differential equation. 

T is the independent' variable t. 

PSI is the value of nth derivative 4> ^ (t) = i|/ft,<j>,<f>», . . 
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PHI is an array of n locations used to store the values of <1> ^ (t) 
(i = 0, 1, . . . , n - 1) in the order of ascending i. 


Other Subroutines Required 

1. AL CHBY 

2. AL DPIN 


Evaluation of Solution and Derivatives 


The finite series 


M+n-i 

B k l)T k (t) Ci ~ 0,1, . . . ,n) 
k=o 

may be evaluated as a function of the independent variable t by means of 
subroutine AL CHBY (see below), which in this case may be accessed via the 
statement 

CALL CHBY (B(l, I + 1) , M + 1 + N - I, T, SUM) 

where 

I, N, M are integers that take on the values of i, n, and M, respectively. 
B(1,I+1) is the location of Bq^ . 

T is the independent variable. 

SUM is the evaluated result. 


Subroutine AL CHBY 

AL CHBY Evaluation of A Finite Series of Chebyshev Polynomials 
FORTRAN IV, Double-Precision Subroutine 

Purpose 

This subroutine is used to evaluate a finite series of the form 

N 

SCt) = ]T)' B k T k (t) , -1 < t < 1 

k=0 
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where (t) are Chebyshev polynomials def ined by 

) , 

T, (t) = cos (k cos -1 t) 

Usage 

This subroutine is entered via the statement 
CALL CHBY (B, N + 1, T, SUM) 

where 

B is a double-precision array of N + 1 locations used to store B^ in the 
order of ascending k. 

N is the order of the highest order Chebyshev polynomial. 

T is the independent variable t and is a double-precision variable. 

SUM is the sum S(t) and a double-precision variable. 

Method (See eq. (A23).) 


Subroutine AL DPIN 


Identification 


AL DPIN, Matrix Inversion 

FORTRAN IV, Double-Precision Subroutine 

Purpose 

This subroutine is used to calculate the inverse of a square matrix A. 
Usage 


This subroutine is entered via the statement 
CALL DPIN (A, N, KDET) 

The parameters are defined as follows: 

A is a double-precision two-dimensional array of dimension N x N used to 
store elements of the matrix A. Upon return, the inverse A -1 will be 
found in the array A. 

N is the order of matrix A. 
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KDET is an error code. KDET = 0 is a normal return; KDET = 1 indicates that 
A is singular. 


Method 

Jordan's method of elimination is used to calculate A -1 (see ref. 14). 
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